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Abstract 

A free-boundary Tokamak Equilibrium Solver (TES), developed for advanced study of tokamak 
equilibra, is described with two distinctive features. One is a generalized method to resolve the 
intrinsic axisymmetric instability, which is encountered after all in equilibrium calculation with a 
free-boundary condition. The other is an extension to deal with a new divertor geometry such 
as snowflake or X divertors. For validations, the uniqueness of a solution is confirmed by the 
independence on variations of computational domain, the mathematical correctness and accuracy 
of equilibrium profiles are checked by a direct comparison with an analytic equilibrium known as a 
generalized Solov’ev equilibrium, and the governing force balance relation is tested by examining 
the intrinsic axisymmetric instabilities. As a valuable application, a snowflake equilibrium that 
requires a second order zero of the poloidal magnetic field is discussed in the circumstance of 
KSTAR coil system. 
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I. INTRODUCTION 


In tokamak physics, plasma equilibrium is a fundamental and essential element to un¬ 
derstand not only the basic equilibrium properties but also various plasma phenomena such 
as MHD instabilities, plasma transport and turbulence, plasma flows and waves, and so on. 
Therefore, various numerical or analytical equilibrium studies [1] have been conducted for a 
long time since the axisymmetric plasma equilibrium relation was established in a general 
form, known as Grad-Shafranov equation m- 

Depending on the characteristics of applications, the studies can be categorized into two 
types of problems. One, so called ’fixed boundary equilibrium’, is solving an equilibrium 
assuming that the plasma boundary or plasma region is known. So, the external equilibrium 
field is ignored and the internal equilibrium profiles and flux distributions are mainly con¬ 
cerned. The other one, so called ’free boundary equilibrium’, is solving the equilibrium with 
unknown plasma boundary. Hence the plasma position and shape (i.e. plasma region) need 
to be obtained as a solution, in addition to those equilibrium profiles and flux distributions. 

Due to the importance of equilibrium as a basis for various physics studies, the majority of 
equilibrium studies has been devoted to the fixed boundary equilibrium problems, while less 
interests given to the free boundary equilibrium solutions. However, recently new demands 
for the free boundary equilibrium analysis have been arisen and turned out to be important. 
For instance, another type of equilibria with new topological features, so called snow-flake 
(SF) divertor [0] or X-divertor [5j equilibria, have been proposed and actively studied in 
various devices l n\ recently. Specially, since the SF divertor configuration requires a 
second-order zero of the poloidal magnetic field, it is now an important issue that should be 
addressed in terms of a free boundary tokamak equilibrium [8]. 

Accordingly, a free boundary Tokamak Equilibrium Solver, called as TES, has been de¬ 
veloped with an emphasis on applications to a design work of plasma equilibrium control 
and to advanced equilibrium study. The developed TES code is featured by two distinctive 
functionalities; a generalized method for stabilization of axisymmetric instabilities, and an 
extension to deal with a second-order zero of the poloidal magnetic field. 

In section [TTJ the numerical solution methods and procedures used in TES is described 
for two types of free boundary equilibrium problems, i.e., ideally free and semi free bound¬ 
ary problems that will be defined therein. Most of numerical techniques and issues have 
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been well known, so that only a brief description for each issues is given unless necessary. 
For validations of TES code, a direct comparison with a generalized analytic solution is 
described in section E3 in addition to the uniqueness of a solution. In section [IVj an intrin¬ 
sic axisymmetric instability, encountered after all during the numerical procedure, is tested 
by examining the variations of plasma equilibria, and a generalized stabilization method is 
introduced and tested. In section [V| an extended feature to deal with a snowflake divertor 


is explained and discussed, followed by a summary and conclusion in section VI 


II. A SOLUTION METHOD IN TES 

A basic numerical method and procedure of an axisymmetric tokamak plasma equilibrium 
with a free boundary condition has been well established [MU- The mathematical and 
numerical treatments used in TES code is also basically in line with those in the references, 
except some improved and extended features. Therefore, in this section, we describes the 
basic numerical treatments and procedures used in TES code briefly unless necessary. 


A. Force balance relation for free boundary plasma equilibrium 


In a toroidally axisymmetric system like a tokamak, the force balance relation of a plasma, 
i.e. plasma equilibrium, can be expressed by a second-order partial differential equation, 
known as Grad-Shafranov equation [21E], using a cylindrical coordinate system ( R , </>, Z) 
with an ignorable (due to axisymmetry) toroidal angle coordinate </>. 


A *ip(R,Z) = -/ 'jl q RJ^ p i(R,Z ) 

nv’Wy 


J<t>,pi{Ri Z) = Rp'(^) + 


ho R 


( 1 ) 

( 2 ) 


where the poloidal flux function (equal to the actual poloidal magnetic flux divided by 2i r) 

is defined by ip(R, Z) = RA^R, Z) from B = V xA (i.e. V • B = 0), the Shafranov operator 

defined by A* = f? 2 V • —, and the prime denotes g' = ——. And p(ib) is an isotropic plasma 

R 2 dip 

pressure, = RB, P a toroidal field function, and J, Pp> i(R. Z) a toroidal current density 

of plasma. Since the J f p pp i{R,, Z) as a source term in Eq. {[Tj) has a strong dependency on 
HR,Z) byEq. (§, it gives rise to a strong non-linearity on the equations. 

In order to deal with a free boundary condition in equilibrium calculation, the Eq. © is 
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generalized by including arbitrary toroidal conductor currents as follows. 


A *ip(R,Z) = ~n 0 RJ^(R, Z) (3) 

J<fr(R,Z) = J ( f >jP i(R,Z) + J^ ,cond(Ri Z) 


where J^ tC emd(R, Z) is the toroidal current density for a conductor. The toroidal conductor 
could be any toroidal current source that can affect the equilibrium force balance, such as 
poloidal field (PF) coil currents or axisymmetric eddy currents on surrounding conductor 
structures. Assuming discrete conductors with uniform current distributions inside, the 
toroidal conductor current density can be expressed by 


J, 


4>,cond 
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k =1 
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( 4 ) 


where J con d,k , hond,k , S k , and fl CO nd,k are the toroidal current density, the toroidal current, 
the cross-sectional area, and the domain region of k-th conductor, respectively. 

Meanwhile, the toroidal current density of plasma J 4pp i in Eq. ([2]) can be set into a 
canonical form ra as shown below 


J<p, P i(R, z) — 
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( 5 ) 


with jfyi'ipaj'ipb) = (1 — , where R geo is the major radius as a reference length scale, 

by the flux per radian at the plasma magnetic axis, py the flux per radian at the plasma 
boundary, and j a suitable profile function. The A and /3 0 are adjustable variables to satisfy 
equilibrium constraints which will be discussed later, while the a m and a n are input variables 
specified by users. Note that J^p^R, Z) is automatically set to zero at the plasma boundary 
in this form by using a normalized poloidal flux, ip s = (-0 — — Vb)- 

In short, a free boundary plasma equilibrium can be obtained by solving Eq. Q with a 
toroidal current density specified by Eqs. @ and (j5j), and corresponding equilibrium profiles 
such as p('tjj) and F('tjj) can be obtained from Eqs. ([2]) and (|5j). 
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B. Numerical approximation by discretizations 


The governing equation described above, i.e. Eq. 0. can be thought as a 2D Poisson’s 
equation in toroidal geometry, so that easily solved using various numerical methods if 
the source term is known. For numerical treatments, the equation is converted to a linear 
algebraic equation by using the centered finite difference method (FDM) [13J on a rectangular 
computational domain in (. R , Z) space, where the grids, (7?/, Zj), are built by 

Rl = Rmin + A.R X (7 — 1), A R = ( Rmax ~ Rmin )/ (NR ~ 1) 

Zj = Z m i n + A Z x (j j — 1), A Z = ( Z max — Z m in)l ( N z — 1) (6) 


with l — 1, • • • , N r and j = 1, • • • , N z . Then the algebraic equation converted by FDM can 
be expressed as follows. 


(AZ) 2 rj 


ipj-1,1 + 
+ 


+ 


(A R) 2 2Ri(AR) 

( 1 




1 + 


+ 


l 


(A R) 2 (AZ) 


tU 


1 


(Vb+D — : ~l- l ()R-iRj>j,i 


(7) 


V(A R) 2 2Ri(AR)J (AZ) 2 

where ip j:i = ip(Ri, Zj) and = J<p(Ri, Zj) with / = 2, • • • ,N R - 1 and j = 2, • • • , N z - 1. 
This algebraic equation can be solved by either using a matrix inversion after reforming 
it in a form of Ax = b or using an iterative method such as multi-grid method [T3J or 
double cyclic reduction na, with an appropriate boundary condition. In TES code, the 
successive-over-relaxation (SOR) method |13j is used as a basic numerical scheme for the 
simplicity. 


C. Iterative solution for non-linearity 


To solve the Eq. Q in the given form, the source term on the right hand side should 
be known. However, the plasma part of the source term has a strong nonlinear dependency 
on ip]d according to Eq. ([2]) or ([5]). To deal with this non-linearity, an iteractive method, 
known as Picard iteration [13], is adopted. Then, the Eq. 0 is expressed as follows 
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where (n) indicates the n-th Picard iteration. Note that the source term in the n-th iteration, 
J^j\, is expressed as a function of i.e. the poloidal flux in the (n-l)th iteration. Hence, 

the is obtained from Eq. (8) using J^ that was evaluated from Then, the J 

is updated using the refreshed ip (n> and provided as a new source term into Eq. (|8]). This 
recursive iteration is continued until a convergence criterion, ||?/4 n ) — ?/4 n_1 ) || < e, is satisfied. 


D. Boundary conditions 


In general, the boundary condition in free boundary equilibrium calculation is not con¬ 
stant and varied due to changes of plasma boundary and equilibrium profiles during the 
numerical iterations, while in a fixed boundary equilibrium it is fixed to zero (ipbc = 0) 
usually. The Dirichlet boundary condition on the edge of a computational domain can be 
provided directly by using a Green’s function formulation [16j. 


<KR,z) 


G(R, Z- R', Z')dR!dZ' 


(9) 


where G(R, Z; R ', Z') is the free space Green’s function which gives the poloidal flux at 
(R, Z) from a unit toroidal current source at ( R',Z '). The free space Green’s function is 
defined by 

G(R, Z,R', Z') = !^y^W[(2-k 2 )K(k)-2E(k)\ 

2 _ 4 RRf 

= (R + i?') 2 + (Z - ~Z/f 

where K(k) and E(k) are elliptic integrals of the first and the second kind [13], respectively. 
Using this, the poloidal flux at the boundary of computational domain can be directly 
obtained by taking into account both plasma and conductor currents as follows 



dlqBi.Z*) = [ G(R b ,Z b] R',Z')J ^( R '’ Z ') d(l ' r i 

Ju>, 

pi 

+ [ G(R b ,Z t -R',Z’)ji:L nd (R',Z')'^' am j (11) 

JQ' , 

cond 

where (Rb, Zb) is the boundary point of the computational domain. Note that Z') 

is varied in every steps of Picard iterations, while Z r ) is not changed unless the 

plasma boundary is specified, which will be discussed later. 
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E. Determination of plasma boundary 


For a stable convergence of the solution, it is important to accurately determine the 
plasma region or boundary in terms of ?/y in every steps of Picard iteration. Generally a 
plasma boundary is formed either by limiters (a limited plasma) or by magnetic fields with 
an X-point (a diverted plasma). Assuming I p > 0, the poloidal flux, ip(R, Z), has a convex 
distribution inside plasma, thus ip a > by. Therefore, the by is defined by the maximum 
value among all poloidal fluxes from limiters and from X-points. When I p < 0, the by is 
defined by the minimum value in a same logic. 

More precisely, both magnetic axis and X-point have a null-field (| V-01 2 = 0), while they 
have different signs of second-derivatives |9j, defined by 

_ / d 2 ip\ f d 2 ip\ / d 2 ip \ 

’ ’ ~ \dR 2 ) ydZ* 2 ) \dRdz) 

If S > 0, the field-null point is a magnetic axis (ip a ). Otherwise (S < 0), it is an X-point. 
The accurate location of the magnetic axis or the X-point is determined by using the Powell’s 
conjugate direction method [13] based on a 2D bicubic interpolation. 

F. Constraints on plasma equilibrium 

In order to have a unique equilibrium solution for Eq. (|2]), a few constraints on plasma 
equilibrium quantaties are necessary. Considering the functional form of Eq. ([5]), two 
constraints, total plasma currents and poloidal plasma beta, are applied. Note that the 
equilibrium constraints could be different when a different functional form of J^ p i(R, Z) is 
used instead of Eq.([5]). For instance, if q(ip) profile is used in J^ p i(R, Z), then q a = (/(by) 
could be used as another appropriate constraint dt 

The constraints can be expressed as 

Ip = j Po)dfl (13a) 

^pi 

a _ Wo)) 

Pp (B 2 UJ 2^0 

where /i 0 = 47 r x 10 _7 [A^/A 2 ] is the permeability of vacuum, a the minor radius, and B p the 
poloidal magnetic field. The braket (•) means an average over a magnetic surface. Therefore, 
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by combining these two equations, the /?o and A can be determined thus giving a unique 
solution. 


G. External equilibrium fields with specified plasma boundary 


In principle, for a free boundary equilibrium problem, the plasma boundary is solved as 
a part of solutions under given external equilibrium fields. In practice, however, it is more 
useful and convenient to solve the equilibrium with a specified plasma boundary. In this 
study, we distinguish them by calling the former as an ideally free boundary problem while 
the latter by a semi free boundary problem. In the case of semi-free boundary, the external 
coil currents are adjusted to provide a required equilibrium field. If a plasma boundary 
is specified in a series of points, the required external equilibrium field currents can be 
determined by solving a minimization problem as shown below 

^'bridry 

( 

^Xpt (^Vcoii 

+ l {pB R (Rj, Zj] Ri , Zi) • A/ coi i^ — B r (Rj , Zj) 

3 =1 l 1=1 

Nxpt ( A^coil 

+ yy l (pB z {Rj, Z^ Ri, Zi) • A/coii,^ — B z (Rj, Zj) 

3=1 l i= 1 



mm 

A/coil 


X 

j = l 



IVcoil 

+ 7 2 ^-fcoil,) 
2=1 


(14) 


where ( Rj,Zj ) is the specified j-th boundary point, A xp(Rj,Zj) = by — tp(Rj,Zj) is the 
poloidal flux error on the point, B R (Rj,Zj ) and B z (R n Zf) are the radial and vertical 
magnetic fields there, and 7 is a Tikhonov parameter for regularization [18]. If an X-point 
is specified as a part of plasma boundary, then the radial and vertical magnetic fields at the 


point should be zeros. This constraint is added as the second and third terms in Eq.(14) 

1 OG 1 3G 

with Gb b = — — 77 = and Gb 7 = + — from this, the external equilibrium held currents 

D A r 7 z, R 

r( n_1 ) 1 A t .. . r( n_1 ) 


RdZ 

r(«) 


are obtained by Ijj ■ = ijj ■ + AI co nj, where Ijj ■ is the coil currents in (n — l)th Picard 
iteration. 
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III. VALIDATIONS OF TES 


According to the numerical methods and procedures described above, a free boundary 
tokamak equilibrium solver (TES) has been developed. For the validations of this code, the 
uniqueness of a solution is firstly checked by examining the independence on the variations 
of computation domains, and the mathematical correctness and accuracy of equilibrium 
profiles are assessed by a direct comparison with an analytic equilibrium solution. 




FIG. 1: (Color online) Two free boundary equilibria, obtained by TES with identical equilibrium 
constraints, are directly compared. On the left, the poloidal magnetic flux in a large computational 
domain (black solid line) is compared with that in a small computational domain (cyon dotted line). 
On the right, several equilibrium profiles are directly compared, such as the isotropic pressure, the 
toroidal field function, the toroidal current density, and the safety factor profile. 


A. Uniqueness of a solution under numerical variations 

Since we are solving the problem in a numerical approach, one fundamental test, which 
is seldom seen in the related literatures, is to examine if it provides an identical result, 
independent on the number of grids or the change of computational domain. Particularly it 
is essential and critical when a free boundary condition is imposed. 

A comparison of two free boundary equilibrium solutions, one in a large and the other 



















in a small computational domains, is shown in Fig. [l] where I p = —2.0 MA, Br = —2.7 T, 
a = 0.48 m, and f3 p = 0.5 with a large elongation k = 2.0. The poloidal magnetic fluxes are 
compared on the left and several equilibrium profiles on the right. The solution for a large 
computational domain (black solid line) was obtained in 0.7 < R < 2.8 m, —1.9 < Z < +1.9 
m with Nr x A+ = 65 x 85, while the one for a small computational domain (cyon dotted 
line) in 1.1 < R < 2.4 m, —1.3 < Z < +1.3 m with Nr x Nz = 45 x 65. As expected, 
the poloidal magnetic fluxes and the equilibrium profiles are shown to be almost identical 
for both. Therefore it confirms that the equilibrium obtained by TES provides a unique 
solution, independently on any change of computational domain and the grid size. 



FIG. 2: (Color online) The poloidal magnetic flux obtained from TES (black solid line) is compared 
directly with the analytic one from the Solovev’s equilibirum (yellow dotted line). The filled contour 
plot with a rectangular boundary shows the full distribution of poloidal magnetic flux including 
vacuum region, which is obtained from TES. 


B. Benchmark with an analytic solution 

For a direct validation of TES, an analytic fixed boundary equilibrium solution, known as 
a generalized Solov’ev equilibrium [l]5j, is considered and compared with a TES result. Note 
that it is to check the mathematical correctness and accuracy of the solution from TES. The 
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pressure and toroidal field function in the analytic solution are assumed to be constant 


~d o 


dp 

dip 


A, 



a 2 


Then, the equilibrium solution can be expressed explicitly as follows 

ip(R, Z) = a + c 2 R 2 + c 3 (.R 4 - 4 R 2 Z 2 ) + c 4 [R 2 In (R) - Z 2 ] + ^-A x - A 2 

o 2 


( 15 ) 


(16) 


where four constants, q, i — 1, • • • ,4 are determined to satisfy the boundary conditions from 
specified plasma boundary, and other two parameters, A\ and A 2 , are adjusted to meet the 
equilibrium constraints. Four boundary conditions, with a modification for the comparison, 
are given by ip(R in , Z in ) = ip(R out , Z out ) = ip(R top , Z top ) = ip b and z } = 0, where 

in, out , and top are the inner-, the outer-, and the top-most boundary points respectively. 
The equilibrium constraints are the total plasma currents I p and the poloidal beta j3 p as 
follows 




4 = 


A 


J J^dRdZ = — 

^ pi 

8i r f p dRdZ 
do 


/ 


I 2 

p 


87r 


— | dRdZ 
do 


\ 




A\ + 


7 


(ipb — ip)dRdZ ) Ai 


/(7s) 


\ 


dRdZ 


Ao 


J 


(17a) 

(17b) 


,W r ? 

thus determining appropriate values of A\ and A 2 . 

The poloidal magnetic fluxes obtained from TES and the analytic solution are directly 
compared in Fig. [2j where I p = 0.5 MA, Bt = 2.7 T, f3 p = 0.5 with elongation k = 1.45 
and minor radius a = 0.5 m in a limited configuration (i.e. without null point). The full 
distribution (including vacuum region) of poloidal magnetic flux from TES is shown as a 
filled contour plot with a rectangular boundary, to show it is indeed a free boundary solution. 
The poloidal magnetic fluxes in plasma region are directly compared by overlapping them; 
one is from TES (black solid line) and the other from the analytic solution (yellow dotted 
line). As shown, two results are not distinguishable and thus the difference is negligible. 
It confirms that the equilibrium informations inside plasma region, obtained by TES, are 
accurately consistent with those from analytic calculations. 

Summarizing two validation results above, it is confirmed that TES provides a unique 
equilibrium solution with high accuracy, consistent with theoretical analysis. 
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IV. AXISYMMETRIC INSTABILITY AND ITS STABILIZATION 


A. Axisymmetric instability of shaped plasma equilibrium 


A tokamak plasma equilibrium has 2D axisymmetric, instrinsic instabilities associated 
with plasma shaping. The most important 2D axisymmetric, i.e. the toroidal mode number 
n=0, instability is known as a vertical instability which becomes unstable once a plasma 
elongation is increased above a threshold. It has been well understood that this instability 
is originated from a J^ p i x B extpo [ force on plasma by an external equilibrium field due to a 
bad curvature associated with the plasma shape. The field curvature can be evaluated by a 
field decay index, ndecay, defined as 


^decay(-fi') 


R dB z 
~B~ z ~dR 


R dB R 
B z dZ 


(18) 


Theoretically, it is well known that a vertically elongated plasma can be unstable when 
^decay < 0 and a radially elongated plasma unstable when nd eC ay >3/2 [20] . 





FIG. 3: (Color online) Decay index of external equilibrium fields (nd eC ay) vs plasma elongation (k) 
for two different plasma betas (/3 P ). Theoretically stable and unstable regimes in terms of nd eC ay 
are marked with filled colors. Additionally a radially elongated plasma with n = 0.7 (left) and a 
vertically elongated plasma with k = 2.0 (right) are shown for ft p = 0.8. 

The relation between plasma shape and field decay index can be seen in FIG. [3] for a 
plasma of I p = 0.5 MA and Br = 2.5 T, where theoretically stable and unstable regimes 
are marked with filled colors. Comparison of two plasmas with different beta ((3 P ) shows 
that the plasma with higher f3 p is less unstable and has wider range of stable k, consistently 
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with theory. Note that the series of equilibria in this figure is obtained by specifying the 
plasma boundary (i.e. as a semi free boundary problem), in order to avoid the axisymmetric 
instability due to the bad curvature. 


B. A generalized stabilization for axisymmetric instabilities 


Due to the axisymmetric instability, the direct solution of plasma equilibrium under given 
external equilibrium field (i.e. as an ideally free boundary problem) has a convergence issue. 
That is, a small deviation of plasma from an equilibrium position is inevitable during a 
numerical iteration, so that the plasma could be drifted and eventually diverged either ra¬ 
dially (when ndecay > 1-5) or vertically (when ridecay < 0.0). In the literature, a conventional 
method to resolve the vertical instability of elongated plasma is simply inserting a feedback 
loop [S] by adding artificial feedback coils which are typically a pair of up-down symmetric 
coils to produce a horizontal magnetic field. In this method, the feedback coil currents are 
adjusted to control the vertical position of magnetic axis to a pre-selected target position 
according to the relation below 

feedback = -sign(Z coil ) {C^ZM - Z torw ,) + C 2 (Z<,% - Z££>)} /, (19) 


where Z target is the desired vertical position of the magnetic axis and Z coi i is the vertical 
position of the control coil. A critical drawback of this method is that the desired vertical 
position of the magnetic axis should be known, prior to obtaining it as a solution from the 
equilibrium calculation. Also the constants Cf and Cf are chosen by trial and error. 


In TES code, the Eq. (19) is modified for a general treatment. Instead of controlling the 


vertical position in a feedback manner, we are eliminating the source of vertical instability 
by compensating Bp field at the center of plasma currents in each steps as following. 


-^feedback 9z 
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R ,vacuum 


D* 

^ R, feedback 


( 20 ) 
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■cur^cur 


where the minus sign indicates a compensation, g z is an adjustable constant, and Bp z 


and B 


R, feedback 


are the radial magnetic fields by external equilibrium conductor currents 


(i.e. vacuum field) and by unit currents of vertical stabilizing coils, respectively. Also note 
that there is no Ip dependency in this method. If g z = 1.0, the exactly same Bp field 
is compensated by the feedback currents. In TES, practically 2.0 < g z < 2.5 is used to 
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ensure a general stabilization by using a up-down symmetric pair of coils, which is set to be 
located radially in the middle of and vertically just outside computational domain. For the 
evaluation of Bn, we use the effective current center (R curr , Z curr ) instead of the magnetic 
axis for a better description of axisymmetric plasma motion as following 


-^cur — Y j Z)dQ p i 

l r 

Zcur J B Jfp pi (B., Z)dVL p i 


( 21 ) 


Note that a generalized method for the radial stabilization is not described here (due to lack 
of practical interest) but also possible in a similar way. 


C. Validation of vertical stability and its stabilization 




FIG. 4: (Color online) Comparisons of plasma displacement responses to the initial perturbations 
of A Z = —0.01 m and A R = —0.01 m for the stability test. Two equilibria with k = 1.0 and 
k = 1.2 are tested for the vertical stability (left) and another two equilibria with n = 1.0 and 
k = 0.8 are tested for the radial stability (right). 

For a validation of the generalized stabilization method described above, we first test the 
validity of force-balance relation solved in TES, by considering the natural axisymmetric 
instability. From FIG. [3j it is obvious that the equilibrium with k = 1.2 is expected to be 
vertically unstable (ri<iecay < 0); while the equilibrium with k = 1.0 to be stable («deca y > 0), 
if the force-balance relation in TES is correct. Similarly, the equilibrium with k = 0.8 is 
expected to be radially unstable (r/decay > 1-5), while the equilibrium with k = 1.0 to be 
stable (ridecay < 1-5). Remind that these equilibria were obtained by specifying the plasma 
boundary and thus produced the required external equilibrium fields as a result (i.e. as a 
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semi free boundary problem). To test the natural vertical stability without any additional 
stabilization, the equilibrium analysis is re-performed as an ideally free boundary problem, 
i.e. by specifying the external equilibrium coil currents which were obtained from the FIG. 
[3| As a seed for vertical or radial instability, a small perturbation is added into the initial 
position of plasma boundary, which is used in 0-th Picard iteration. 

The comparisons of vertical and radial displacement responses to small deviations of 
A Z = —0.01 m and A R = —0.01 m are shown in FIG. |4j On the left, the vertical stability 
is tested by an initial perturbation, A Z = —0.01 m, for two equilibria; one with k = 1.0 
(black dotted line) and the other with k = 1.2 (blue dotted line). Consistently with the 
theoretical expectations, the initial perturbation of the former was naturally stabilized, while 
the one of the latter was exponentially diverged. On the right, the radial stability is tested 
by an initial perturbation, A R = —0.01 m, for two equilibria; one with n = 1.0 (black dotted 
line) and the other with k = 0.8 (blue dotted line). Similarly, the initial perturbation of 
the former was naturally stabilized or stable, while the one of the latter was exponentially 
diverged. Therefore, it confirms that the force-balance relation used in TES is correctly 
solved and the associated instability is precisely consistent with the theory. 



FIG. 5: (Color online) Vertical displacement responses to an initial perturbation of A Z = —0.03 m 
are compared for the equilibrium shown on the left. The first one (black) is during TES iteration 
as a semi free boundary problem, the second one (blue) as an ideally free boundary problem, and 
the third one (red) with the generalized vertical stabilization. 
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For the validation of the generalized stabilization method, a further strongly shaped and 
up-down asymmetric plasma is considered as a worst case. The reference equilibrium is 
obtained in a single null (SN) configuration with Ip = —0.5 MA, Bp = 2.5 T, ftp — 0.1, 
and k = 1.7 («decay = — 1.8) as shown in FIG. [5| The comparison of vertical displacement 
responses with and without the generalized stabilization is shown on the right of the FIG. 
[5j The evolution of a semi free boundary solution (black line) shows that it converged 
to Z mag = —0.03 m (thus it is a reference equilibrium position). In case of ideally free 
boundary solution (blue line) without any stabilization, it was slowly drifted upward and 
finally diverged, as expected. Then, by applying the generalized vertical stabilization (red 
line), the evolution was really stabilized so that it was smoothly evolved and converged to 
the reference position closely. Here, the final difference of vertical position compared with 
the reference is about 0.5 cm. Therefore, it demonstrates that the generalized method can 
effectively stabilize the natural vertical instability of elongated plasmas and automatically 
guide the plasma to an equilibrium position, that is consistent with that from a semi free 
boundary solution. 


V. EXTENSION TO ADVANCED EQUILIBRIUM ANALYSIS 

Recently new types of tokamak equilibria have been proposed and studied in various 
devices, in order to resolve the issue of an excessive heat and particle fluxes onto the plasma 
facing components in ITER and beyond. These are featured by a new divertor configuration 
such as snowflake [1] and (super) X divertors [5j. Particularly the snowflake equilibrium 
requires to have a second-order zero of poloidal flux at the null-field point so that it is not 
straight-forward to deal with it by a conventional free boundary equilibrium solver |8J. To 
solve this new equilibrium with specified plasma boundary, the minimization constraint, Eq. 


(14), for required external equilibrium field currents is modified in TES as follows 
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snowflake equilibrium that requires a second-order zero of ^(R, Z) can be directly obtained 
without any special treatment in TES. 




FIG. 6: (Color online) A comparison of two distinctive plasma equilibria. One is with a typical 
double null divertor (left) and the other with a snowflake divertor (right). The plasma boundary 
points specified in calculation are marked with a red circle. 

Figure [6] shows a comparison of two equilibria obtained by TES with identical plasma 
equilibrium parameters, which are Ip =1.0 MA, Bp =2.5 T, k =2.0, and f3p =0.2. One (on 
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TABLE I: Required external coil currents (kA) to form the equilibria shown in FIG. [6] 



Double-Null (DN) 

Snow-Flake (SF) 

PF1 

-6.02 

-1.04 

PF2 

9.64 

-19.05 

PF3 

6.10 

135.27 

PF4 

7.07 

-67.32 

PF5 

8.33 

18.60 

PF6 

-1.03 

-5.30 

PF7 

-7.26 

-4.03 


the left) is a typical double null (DN) divertor and the other (on the right) is a snowflake 
(SF) divertor configurations. The difference of two equilibria is easily seen from the magnetic 
distributions around the field null points. In the SF configuration, it is clearly seen that 
three concave and another three convex distributions are formed alternately, centred at the 
up-down symmetric field null points, i.e. a second order zero of poloidal magnetic flux is 
formed. 

Table [T] shows the external equilibrium coil currents required to form the target equilibria 
shown in FIG. [6| It is important to note that in the case of SF equilibrium some of coil 
currents are required extremely large values, while in the case of DN equilibrium all coil 
currents are well balanced. It indicates that it is not practically possible to form the SF 
equilibria onto the KSTAR by using current coil system. Therefore, a new coil system, 
specially designed for SF divertor, is essentially needed. In fact, it is consistent with the 
recent highlighted issue |8] in the study on advanced divertor configurations. In addition, it 
is worthwhile to note that the SF equilibrium here is solved self-consistently by considering 
full force-balance relations in a toroidal system, while in the reference [8], it is solved by 
using a simplified wire plasma model. 

VI. CONCLUSIONS 

A free-boundary tokamak equilibrium solver, developed for advanced study of tokamak 
equilibra, was described with various validation results. The developed solver, named as 
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TES, is characterized by two distinctive features. At first, a generalized stabilization method 
for intrinsic axisymmetric instabilities was applied, which is encountered after all in equilib¬ 
rium calculation under a free boundary condition. In this method, the source of axisymmet¬ 
ric instabilities is directly removed or minimized, instead of feedback controlling the plasma 
position to a target location. Thus, it ensures in general that the TES code produces a 
solution stably even under highly (axisymmetrically) unstable conditions. 

The other important feature is an extension to deal with a new divertor geometry such 
as snowflake or X divertors. To deal with the innovative divertor concept, particularly the 
snowflake divertor, the equilibrium solver needs to be able to control the location of second 
order zero of poloidal magnetic field. By implementing this functionality into the TES 
code, it was demonstrated that the snowflake type of advanced tokamak equilibria can be 
analysed in consideration of full toroidal force balance relations, instead of using a simplified 
wire plasma model. 

For the validation of TES code, the uniqueness of a solution was confirmed by the inde¬ 
pendence on variations of computational domain, the mathematical correctness and accuracy 
of equilibrium profiles were checked by a direct comparison with the generalized Solov’ev 
equilibrium, and the governing force balance relation was tested by examining the intrinsice 
axisymmetric instabilities. 

As a valuable application, a snowflake equilibrium was analysed by taking into account 
the KSTAR equilibrium coil system. Since the KSTAR has a limited set of equilibrium 
control coils, it is important to check whether the innovative divertor equilibria can be 
realized in the current system. The analysis results suggest that practically it is not possible 
to form a snowflake equilibrium in current KSTAR device so that additional control coils 
need to be considered for the study of advanced divertors in future. 
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